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Abstract 

Nonlinear coupling between inter- and intra-element dynamics appears as a collective behaviour 
of elements. The elements in this paper denote symptoms such as a bacterium having an internal 
network of genes and proteins, a reactive droplet, a neuron in networks, etc. In order to elucidate 
the capability of such systems, a simple and reasonable model is derived. This model exhibits the 
rich patterns of systems such as cell membrane, cell fusion, cell growing, cell division, firework, 
branch, and clustered clusters (self-organized hierarchical structure, modular network). This model 
is extremely simple yet powerful; therefore, it is expected to impact several disciplines. 
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We consider non-equilibrium elements whose internal dynamics interacts with their 
macroscopic or mesoscopic order(i-5). An assembly of such elements can describe one aspect 
of soft-matter; there are exciting projects concerning soft matter in several countries includ- 
ing Japan, where a huge project started recently. These reveal the world-wide importance 
of such studies. In this paper, the elements denote symptoms such as a bacterium having an 
internal network of genes and proteins(^), a reactive droplet in reaction-diffusion systems(7), 
a neuron in networks (8, 9), etc(10-14)- These elements exhibit not only spatio-temporal pat- 
terns but also collective functions. For instance, the cohort migration of mammalian cells 
forms tissue patterns (1 5,16), and the Proteus mirabilis effectively invades human urothelial 
cells by swarming(i 7). Further, swarm intelligence has been extensively studied in order to 
enable a collection of simple robots to perform advanced tasks (is). Here, we show a simple 
model derived by means of mathematical techniques to study the cross-cutting phenomenon 
underlying the above systems while ignoring system-specific details. The derived model ex- 
hibits rich patterns such as a modular network and a closed membrane that moves around, 
grows, and multiplies like a cell. This model is expected to shed light on issues related to 
several disciplines. 

In order to avoid an ad-hoc toy model, we propose a broad class of models from which 
we derive a simple model by means of centre-manifold reduction and phase reduction(i5). 
At the end of this paper, we briefly discuss the vast possibilities of this derived model. 

Self-sustained (or limit-cycle) oscillator is the leading candidate for the simplest dy- 
namical element. Thus, we assume a supercritical Hopf bifurcation for the intra-element 
dynamics. A simple interaction among the elements is mediated by a chemical that diffuses 
in space. Thus, we assume that the elements exhibit chemotaxis(^f), which means that 
the elements are driven by the gradient of chemical density, and the elements produce and 
consume the chemical depending on their state. We carry out centre manifold reduction 
in the neighbourhood of the Hopf bifurcation point. However, because the reduced model 
captures the critical centre of systems, the model describes the systems in a broad param- 
eter space until another bifurcation occurs. In addition, the reduced model can be derived 
from another class of models having the same critical centre as that in our original class of 
models. Thus, the model derived in the following is fairly universal. 
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The reduced model of the chemotaxis oscillators is given 

fa = l + {KP(ri) + c.c), (1) 
r. = -VP(r)| r=ri + c.c, (2) 

4>i represents the internal state of ith element. The D-dimensional real vector represents 
the position of the ith element, k is a complex constant. P represents an interaction among 
the elements and is given by 

P(r) = ^e-^^^-r), (3) 
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where c is a real parameter. Note that the model becomes a potential system when ip is 
adiabatically eliminated. In one-dimensional space, the coupling kernel G is simply expressed 
as G D=1 (r) = j-e~ p l r l. b and p are complex constants. In two-dimensional space, G D=2 (r) = 
■^K (p\r\), where K is the modified Bessel function of the second kind with a complex 
argument. In an any-dimensional space, G(r) oscillates and rapidly decreases when \r\ 
increases, and G(r) almost vanishes when \r\ is greater than the coupling length r c = 1/Rep. 
Thus, the main characteristic of G(r) - oscillating and decreasing - is qualitatively well 
described by G D=1 (r), which we substitute for G in an any-dimensional space for simplicity. 
In fact, we numerically confirmed that the spatio-temporal patterns shown in this paper can 
be observed for the original G with slight and suitable parameter changes. 

We now rescale Eq. (CQ) and Eq. (j2j). Introducing a variable ipi defined as ipi = c{<pi — 
[1 + (kG(0) + c.c.)]t}, we rescale the spatio-temporal coordinate as r' = Repr« and d t > = 
\p/(cnb)\d t . Omitting the prime, we obtain 

tpi = e ~ lBdil Sill (^ + a \ R 3i\ ~ C 0> ( 4 ) 

n = c 3 ^2 Rjie' 111 ^ sm(V ji + a\Rji\ - c 2 ), (5) 

where Rji = rj — Rji = Rji/\Rji\ and = ipj — ipi. These equations contain the 
four real parameters c\ = a.rg(cnb/p) — ir/2, c% = arg(— b) — tt/2, c 3 = Rep\p/cn\(> 0) and 
a = Imp/Rep(> 0). Note that C3 is the ratio of the time scales of ipi and tv 

Hence, we derived two models: model I Eq. ([1]) and Eq. (j2J), and model II Eq. 01]) and 
Eq. (j3J). (The second model is equivalent to the first model in one-dimensional space. In 
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a higher-dimensional space, the former is an approximation of the latter as stated above.) 
These models are extended models reported in the previous paper (2 1,22). 

Now, we show the richness of these models by using numerical calculations carried out in 
two-dimensional space with a periodic boundary condition. The boundary conditions are not 
important in sufficiently large systems because the coupling function G(r) decays rapidly 
as | t* j increases. The initial condition is such that the positions and phases are randomly 
distributed. The number of elements is fifty. The figures show snapshots of the element 
distribution in two-dimensional space after the transient time. The colours represent the 
phase (or ip) of the element. Figurel shows firework exhibited by the first model. This 
pattern is static in space, and the phase waves spread from the centre of the pattern, which 
can correspond to the target pattern in reaction-diffusion systems. In fact, our models also 
exhibit spiral waves with another choice of parameters. Figure2 shows the closed membrane 
exhibited by the second model. This membrane moves around while maintaining almost 
its shape. We can observe the fraction of the membrane where the elements are relatively 
synchronised ( 23). This fraction is puffs up as time passes. Thus, in this figure, the mem- 
brane subsequently moves upwards. If we change the parameters slightly, the membrane 
is divided into two synchronous clusters, which grow to form two closed membranes, and 
this process repeats. When the number of elements constituting one membrane decreases 
due to membrane division, the elements merge with neighbouring elements. This dynamical 
pattern reminds us of the proliferation of cells. Figure3 shows clustered clusters (or modular 
networks) exhibited by the second model. This self-organised hierarchical structure is consti- 
tuted by synchronous clusters that exhibit anti-phase synchronisation with the neighbouring 
clusters, which reminds us of the self-differentiation of cells. Here, it should be noted that 
in these three patterns, adjacent elements exhibit an approximate in-phase synchronisation. 
We cannot present all the patterns in this paper; however, by simply changing the parame- 
ters, our models can exhibit a junction of three branches, a crystalline lattice, gas, collective 
translational motion parallel to plane phase-wave, stick-slip motion of clusters, train motion, 
etc. 

In order to understand the patterns analytically, we consider two-oscillators system; this 
corresponds to the case where there is only one oscillator in the neighbourhood of the ith 
oscillator. This analysis sheds light on many-oscillators system too. Because ri\\Rji, the 
two oscillators move only along a line parallel to Rji that does not change with time. Thus, 
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FIG. 1: Firework. Snapshot of element distribution in two-dimensional space. The colours repre- 
sent the internal state (ft of the elements. The parameters are k = 0.8 + il, p = 0.1, 6 = 2 + i3.5 
and c = 1. The space size is 120 x 120, and it is shown in entirety. Although we have adopted a 
point element in this paper, we plot its position with a finite size for visualizing. 




FIG. 2: Closed membrane shown in the same manner as Fig.l. The parameters are c\ = 1.5, C2 = 3, 
C3 = 0.02 and a = 0. The space size is 10 x 10. 



we use Ti = T*i • Rji instead of ty The difference of the two oscillators ^ = ip2 ~ ipi and 
R = r2 — r\ obeys \& = — 2e~ R cos(aR — c\) sin R = 2c^e~ R sin(aR — C2) cos where we 
assume R > without losing generality because we can transpose the labels of oscillators 1 
and 2. When the first equation is divided by the second equation, we can separate variables 
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FIG. 3: Clustered clusters (or Modular networks) shown in the same manner as Fig.l. The param- 
eters are c\ = c 2 = C3 = 1.5 and a = 1.6. The space size is 30 x 30. 

* and R; then, integrating once, we derive the invariant curve 

I sin* I = Ee aiR \sm(aR- c 2 )\ a \ (6) 

where a\ = sin(ci — c 2 )/c 3 and a 2 = cos(ci — c 2 )/(ac 3 ). (We can easily derive another curve 
when a = or C3 = 0.) E is a conserved quantity that is defined by the initial conditions 
*(0) and R(0). Thus, the difference of the two oscillators moves on this invariant curve. 
This equation implies that if the two oscillators synchronise (in-phase * = or anti-phase 

* = 7r), the distance must be R = c 2 /a mod (it /at)- In fact, this distance can be observed 
in Fig. 3, where the distance between the neighbouring elements in the synchronous cluster 
is c 2 /a, and the distance between neighbouring synchronous clusters is c 2 /a + 7r/a. Further, 
Eq. prohibits the distance from making the right hand side of Eq. become greater than 
one. This implies an effective excluded volume, i.e. the elements spontaneously maintain a 
finite distance between each other even if they do not have an excluded volume. Another 
analysis can also be carried out, and we will show this elsewhere. 

At the end of this paper, we comment on the vast possibilities of our models. (1) These 
models pertain to not only the systems stated in the introductory part but also the following 
(a)-(f). (a) Networks, where the strength of the link (or edge) between nodes corresponds 
to the distance \R%j\ as shown in Fig. Ufa). Actually, closed membrane and clusterd clusters 
correspond self-organised small- world network and modular network respectively. In a social 
networking service, the update of an individual page is as frequent as that of the relation 
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between individual pages, i.e. the time scale of a node is comparable to that of a link. This 
is rarely observed in conventional home pages. Our models can shed light on this type of 
recent network types. In addition, our models suggest new possibilities for information or 
memory processing in neural networks exhibiting spike timing dependent synaptic plastic- 
ity(STDP)(£^) etc. (b) Fluids where an acoustic wave mediates interaction among the radii 
of bubbles, (c) Motile-spin glass, where the state of the spin corresponds to the phase ip 
(or <p), and the wave function of the electron surrounding the spin corresponds to the oscil- 
latory coupling function G{25) as shown in Fig. H^b). (d) Reaction-diffusion systems. For 
instance, the phase waves on the firework in our models correspond to the target and spiral 
patterns in reaction-diffusion systems. Further, the firework, branch, junction of the three 
branches and cell division in our models can correspond to the spots, stripes, defects in the 
stripes and pulse division in reaction-diffusion systems(£#). (e) Frustration systems. In our 
model with some parameters, the elements stabilize when they are located away from each 
other. This causes frustration in high-density elements systems, (f) Time-delayed systems, 
r / implies an effective delay in the interaction among elements. (2) Thermodynamical 
limit, i.e. the presence of many elements, which involves non-equilibrium statistics and an 
extended kinetic theory of gases(^T). We can derive a continuous model for the density of 
elements. (3) The coulomb interaction, i.e. galvanotaxis instead of chemotaxis(^5). This 
will shed light on ionic fluids where dipoles may correspond to the elements. Further, this 
is interesting because of the fact that cancer cells and, more generally, biological cells ex- 
hibit galvanotaxis. (4) Three- or higher-dimensional spaces. We may observe self-organised 
spherical shells. (5) Design of materials having newly identified physical properties by using 
the self-organised structures of our mo dels ( 29,30). 

There are previous studies (referred in the introductory part) sharing similar concepts 
as in this paper. However, to date, we have not come across models that are as simple, 
reasonable, powerful and that provide for analytical results. Such models are derived in this 
paper. Our models show that only one variable (phase in this paper) other than the position 
variable is sufficient to obtain the rich collective behaviour of systems. 
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